module mod_clm_slakehydrology
  !
  use mod_intkinds
  use mod_realkinds
  use mod_runparams, only : dtsrf
  !
  ! Calculation of Lake Hydrology.
  ! Full hydrology, aerosol deposition, etc. of snow layers is done.
  ! However, there is no infiltration, and the water budget is balanced with
  ! qflx_qrgwl.
  ! Lake water mass is kept constant. The soil is simply maintained at
  ! volumetric saturation if ice melting frees up pore space. Likewise, if
  ! the water portion alone at some point exceeds pore capacity, it is
  ! reduced. This is consistent with the possibility of initializing the
  ! soil layer with excess ice.
  !
  ! If snow layers are present over an unfrozen lake, and the top layer of
  ! the lake is capable of absorbing the latent heat without going below
  ! freezing, the snow-water is runoff and the latent heat is subtracted
  ! from the lake.
  !
  ! Minimum snow layer thickness for lakes has been increased to avoid
  ! instabilities with 30 min timestep.
  ! Also frost / dew is prevented from being added to top snow layers that
  ! have already melted during the phase change step.
  !
  implicit none

  private

  save

  public :: SLakeHydrology        ! Calculates soil/snow hydrology

  contains
  !
  ! WARNING: This subroutine assumes lake columns have one and only one pft.
  !
  ! Sequence is:
  !  SLakeHydrology:
  !    Do needed tasks from Hydrology1, Biogeophysics2, & top of Hydrology2.
  !    -> SnowWater:          change of snow mass and snow water onto soil
  !    -> SnowCompaction:     compaction of snow layers
  !    -> CombineSnowLayers:  combine snow layers that are thinner than minimum
  !    -> DivideSnowLayers:   subdivide snow layers that are thicker than
  !                            maximum
  ! Add water to soil if melting has left it with open pore space.
  ! If snow layers are found above a lake with unfrozen top layer, whose top
  ! layer has enough heat to melt all the snow ice without freezing, do so
  ! and eliminate the snow layers.
  ! Cleanup and do water balance.
  ! Do SNICAR stuff and diagnostics.
  ! Call SnowAge_grain (it must be done here because the snow filters at the
  ! driver level are non-lakec only.
  !
  subroutine SLakeHydrology(lbc, ubc, lbp, ubp, num_lakec, filter_lakec, &
                            num_lakep, filter_lakep)
    ! Snow filter for lakes is not returned to driver.  That's okay, because
    ! it looks like it is only needed for the call to SnowAge_grain, which
    ! will be done at the bottom of this module.
    use mod_clm_type
    use mod_clm_atmlnd, only : clm_a2l
    use mod_clm_varcon, only : denh2o, denice, spval, hfus, tfrz, cpliq, cpice
    use mod_clm_slakecon, only : lsadz
    use mod_clm_varpar, only : nlevsno, nlevgrnd, nlevsoi
    use mod_clm_snowhydrology, only : SnowCompaction, CombineSnowLayers, &
                                 SnowWater, BuildSnowFilter
    use mod_clm_snowhydrology, only : DivideSnowLayers_Lake
    use mod_clm_snicar, only : SnowAge_grain, snw_rds_min
    implicit none
    integer(ik4), intent(in) :: lbc, ubc                  ! column bounds
    integer(ik4), intent(in) :: lbp, ubp                  ! pft bounds
    ! number of column lake points in column filter
    integer(ik4), intent(in) :: num_lakec
    ! column filter for lake points
    integer(ik4), intent(in) :: filter_lakec(ubc-lbc+1)
    ! number of pft lake points in column filter
    integer(ik4), intent(in) :: num_lakep
    ! pft filter for lake points
    integer(ik4), intent(in) :: filter_lakep(ubp-lbp+1)

    real(rk8), pointer, contiguous :: frac_sno_eff(:)  ! needed for snicar code
    real(rk8), pointer, contiguous :: qflx_floodg(:) ! gridcell flux of flood water from RTM
    real(rk8), pointer, contiguous :: qflx_floodc(:) ! column flux of flood water from RTM
    real(rk8), pointer, contiguous :: frost_table(:) ! frost table depth (m)
    real(rk8), pointer, contiguous :: zwt_perched(:) ! perched water table depth (m)
    ! perched wt sub-surface runoff (mm H2O /s)
    real(rk8), pointer, contiguous :: qflx_drain_perched(:)
    real(rk8), pointer, contiguous :: qflx_h2osfc_surf(:) ! surface water runoff (mm H2O /s)
    real(rk8), pointer, contiguous :: qflx_snow_melt(:)! net snow melt
    real(rk8), pointer, contiguous :: qflx_rsub_sat(:) !soil saturation excess [mm h2o/s]
    integer(ik4), pointer, contiguous :: pcolumn(:)   ! pft's column index
    integer(ik4), pointer, contiguous :: pgridcell(:) ! pft's gridcell index
    integer(ik4), pointer, contiguous :: cgridcell(:) ! column's gridcell
    integer(ik4), pointer, contiguous :: clandunit(:) ! column's landunit
    ! volumetric soil water at saturation (porosity)
    real(rk8), pointer, contiguous :: watsat(:,:)
    real(rk8), pointer, contiguous :: z(:,:)       ! layer depth  (m)
    real(rk8), pointer, contiguous :: dz_lake(:,:) ! layer thickness for lake (m)
    real(rk8), pointer, contiguous :: forc_rain(:) ! rain rate [mm/s]
    real(rk8), pointer, contiguous :: forc_snow(:) ! snow rate [mm/s]
    real(rk8), pointer, contiguous :: begwb(:)     ! water mass begining of the time step
    ! qflx_evap_soi + qflx_evap_can + qflx_tran_veg
    real(rk8), pointer, contiguous :: qflx_evap_tot(:)
    real(rk8), pointer, contiguous :: forc_t(:)    ! atmospheric temperature (Kelvin)
    logical, pointer, contiguous :: do_capsnow(:) ! true => do snow capping
    real(rk8), pointer, contiguous :: t_grnd(:)    ! ground temperature (Kelvin)
    ! soil evaporation (mm H2O/s) (+ = to atm)
    real(rk8), pointer, contiguous :: qflx_evap_soi(:)

    real(rk8), pointer, contiguous :: dz(:,:)     ! layer thickness depth (m)
    real(rk8), pointer, contiguous :: zi(:,:)     ! interface depth (m)
    integer(ik4), pointer, contiguous :: snl(:)  ! number of snow layers
    real(rk8), pointer, contiguous :: h2osno(:)   ! snow water (mm H2O)
    real(rk8), pointer, contiguous :: snow_depth(:) ! snow height (m)
    ! mass fraction of lake layer that is frozen
    real(rk8), pointer, contiguous :: lake_icefrac(:,:)
    real(rk8), pointer, contiguous :: t_lake(:,:)      ! lake temperature (Kelvin)
    real(rk8), pointer, contiguous :: qflx_snomelt(:)  ! snow melt (mm H2O /s)
    real(rk8), pointer, contiguous :: eflx_snomelt(:)  ! snow melt heat flux (W/m**2)
    ! total sensible heat flux (W/m**2) [+ to atm]
    real(rk8), pointer, contiguous :: eflx_sh_tot(:)
    ! sensible heat flux from ground (W/m**2) [+ to atm]
    real(rk8), pointer, contiguous :: eflx_sh_grnd(:)
    ! heat flux into snow / lake (W/m**2) [+ = into soil]
    real(rk8), pointer, contiguous :: eflx_soil_grnd(:)
    ! Here this includes the whole lake radiation absorbed.
    ! net heat flux into ground (W/m**2)
    real(rk8), pointer, contiguous :: eflx_gnet(:)
    ! net heat flux into lake / snow surface, excluding light
    ! transmission (W/m**2)
    real(rk8), pointer, contiguous :: eflx_grnd_lake(:)

    real(rk8), pointer, contiguous :: endwb(:)         ! water mass end of the time step
    real(rk8), pointer, contiguous :: snowice(:)       ! average snow ice lens
    real(rk8), pointer, contiguous :: snowliq(:)       ! average snow liquid water
    real(rk8), pointer, contiguous :: t_soisno(:,:)    ! snow temperature (Kelvin)
    real(rk8), pointer, contiguous :: h2osoi_ice(:,:)  ! ice lens (kg/m2)
    real(rk8), pointer, contiguous :: h2osoi_liq(:,:)  ! liquid water (kg/m2)
    real(rk8), pointer, contiguous :: h2osoi_vol(:,:)  ! volumetric soil water [m3/m3]
    real(rk8), pointer, contiguous :: qflx_drain(:)    ! sub-surface runoff (mm H2O /s)
    real(rk8), pointer, contiguous :: qflx_surf(:)     ! surface runoff (mm H2O /s)
    real(rk8), pointer, contiguous :: qflx_infl(:)     ! infiltration (mm H2O /s)
    ! qflx_surf at glaciers, wetlands, lakes
    real(rk8), pointer, contiguous :: qflx_qrgwl(:)
    ! total runoff (qflx_drain+qflx_surf+qflx_qrgwl) (mm H2O /s)
    real(rk8), pointer, contiguous :: qflx_runoff(:)
    real(rk8), pointer, contiguous :: qcharge(:)       ! aquifer recharge rate (mm/s)
    ! net water input into soil from top (mm/s)
    real(rk8), pointer, contiguous :: qflx_top_soil(:)
    ! liquid water + ice from layer above soil to top soil layer or
    ! sent to qflx_qrgwl (mm H2O/s)
    real(rk8), pointer, contiguous :: qflx_sl_top_soil(:)
    ! water onto ground including canopy runoff [kg/(m2 s)]
    real(rk8), pointer, contiguous :: qflx_prec_grnd(:)
    ! water onto ground including canopy runoff [kg/(m2 s)]
    real(rk8), pointer, contiguous :: qflx_prec_grnd_col(:)
    ! snow on ground after interception (mm H2O/s) [+]
    real(rk8), pointer, contiguous :: qflx_snow_grnd_pft(:)
    ! snow on ground after interception (mm H2O/s) [+]
    real(rk8), pointer, contiguous :: qflx_snow_grnd_col(:)
    ! rain on ground after interception (mm H2O/s) [+]
    real(rk8), pointer, contiguous :: qflx_rain_grnd(:)
    ! fraction of ice relative to the tot water
    real(rk8), pointer, contiguous :: frac_iceold(:,:)
    ! pft quantity averaged to the column (assuming one pft)
    real(rk8), pointer, contiguous :: qflx_evap_tot_col(:)
    ! factor that reduces ground saturated specific humidity (-)
    real(rk8) ,pointer, contiguous :: soilalpha(:)
    real(rk8), pointer, contiguous :: zwt(:)    ! water table depth
    ! fractional area with water table at surface
    real(rk8), pointer, contiguous :: fcov(:)
    ! fractional area with water table at surface
    real(rk8), pointer, contiguous :: fsat(:)
    ! effective fraction of roots in each soil layer
    real(rk8), pointer, contiguous :: rootr_column(:,:)
    ! ground surface evaporation rate (mm H2O/s) [+]
    real(rk8), pointer, contiguous :: qflx_evap_grnd(:)
    ! ground surface evaporation rate (mm H2O/s) [+]
    real(rk8), pointer, contiguous :: qflx_evap_grnd_col(:)
    ! sublimation rate from snow pack (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_sub_snow(:)
    ! sublimation rate from snow pack (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_sub_snow_col(:)
    ! surface dew added to snow pack (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_dew_snow(:)
    ! surface dew added to snow pack (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_dew_snow_col(:)
    ! ground surface dew formation (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_dew_grnd(:)
    ! ground surface dew formation (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_dew_grnd_col(:)
    ! rain on ground after interception (mm H2O/s) [+]
    real(rk8), pointer, contiguous :: qflx_rain_grnd_col(:)
    ! excess snowfall due to snow capping (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_snwcp_ice_col(:)
    ! excess snowfall due to snow capping (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_snwcp_ice(:)
    ! excess rainfall due to snow capping (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_snwcp_liq_col(:)
    ! excess rainfall due to snow capping (mm H2O /s) [+]
    real(rk8), pointer, contiguous :: qflx_snwcp_liq(:)
    ! irrigation flux (mm H2O /s)
    real(rk8), pointer, contiguous :: qflx_irrig(:)
    !New SNICAR variables from Hydrology1
    ! effective snow grain radius (col,lyr) [microns, m^-6]
    real(rk8), pointer, contiguous :: snw_rds(:,:)
    ! mass of hydrophobic BC in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_bcpho(:,:)
    ! mass of hydrophilic BC in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_bcphi(:,:)
    ! total mass of BC in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_bctot(:,:)
    ! total column mass of BC in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_bc_col(:)
    ! total top-layer mass of BC (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_bc_top(:)
    ! mass of hydrophobic OC in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_ocpho(:,:)
    ! mass of hydrophilic OC in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_ocphi(:,:)
    ! total mass of OC in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_octot(:,:)
    ! total column mass of OC in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_oc_col(:)
    ! total top-layer mass of OC (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_oc_top(:)
    ! mass of dust species 1 in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_dst1(:,:)
    ! mass of dust species 2 in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_dst2(:,:)
    ! mass of dust species 3 in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_dst3(:,:)
    ! mass of dust species 4 in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_dst4(:,:)
    ! total mass of dust in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_dsttot(:,:)
    ! total column mass of dust in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_dst_col(:)
    ! total top-layer mass of dust in snow (col,lyr) [kg]
    real(rk8), pointer, contiguous :: mss_dst_top(:)
    !Additional SNICAR variables from Hydrology2
    ! mass concentration of BC species 1 (col,lyr) [kg/kg]
    real(rk8), pointer, contiguous :: mss_cnc_bcphi(:,:)
    ! mass concentration of BC species 2 (col,lyr) [kg/kg]
    real(rk8), pointer, contiguous :: mss_cnc_bcpho(:,:)
    ! mass concentration of OC species 1 (col,lyr) [kg/kg]
    real(rk8), pointer, contiguous :: mss_cnc_ocphi(:,:)
    ! mass concentration of OC species 2 (col,lyr) [kg/kg]
    real(rk8), pointer, contiguous :: mss_cnc_ocpho(:,:)
    ! mass concentration of dust species 1 (col,lyr) [kg/kg]
    real(rk8), pointer, contiguous :: mss_cnc_dst1(:,:)
    ! mass concentration of dust species 2 (col,lyr) [kg/kg]
    real(rk8), pointer, contiguous :: mss_cnc_dst2(:,:)
    ! mass concentration of dust species 3 (col,lyr) [kg/kg]
    real(rk8), pointer, contiguous :: mss_cnc_dst3(:,:)
    ! mass concentration of dust species 4 (col,lyr) [kg/kg]
    real(rk8), pointer, contiguous :: mss_cnc_dst4(:,:)

    ! New Diagnostics
    ! snow temperature in top layer (col) [K]
    real(rk8), pointer, contiguous :: snot_top(:)
    ! temperature gradient in top layer (col) [K m-1]
    real(rk8), pointer, contiguous :: dTdz_top(:)
    ! effective snow grain size, top layer(col) [microns]
    real(rk8), pointer, contiguous :: snw_rds_top(:)
    ! liquid water fraction in top snow layer (col) [frc]
    real(rk8), pointer, contiguous :: sno_liq_top(:)
    ! mass of snow in top layer (col) [kg]
    real(rk8), pointer, contiguous :: h2osno_top(:)

    integer(ik4)  :: p,fp,g,c,j,fc,jtop   ! indices
    integer(ik4)  :: num_shlakesnowc      ! number of column snow points
    ! column filter for snow points
    integer(ik4)  :: filter_shlakesnowc(ubc-lbc+1)
    ! number of column non-snow points
    integer(ik4)  :: num_shlakenosnowc
    ! column filter for non-snow points
    integer(ik4)  :: filter_shlakenosnowc(ubc-lbc+1)
    ! flag when new snow node is set, (1=yes, 0=no)
    integer(ik4)  :: newnode
    ! layer thickness rate change due to precipitation [mm/s]
    real(rk8) :: dz_snowf
    ! bulk density of newly fallen dry snow [kg/m3]
    real(rk8) :: bifall
    ! snow precipitation incident on ground [mm/s]
    real(rk8) :: qflx_prec_grnd_snow(lbp:ubp)
    ! rain precipitation incident on ground [mm/s]
    real(rk8) :: qflx_prec_grnd_rain(lbp:ubp)
    ! temporary evap_soi limited by top snow layer content [mm/s]
    real(rk8) :: qflx_evap_soi_lim
    ! temporary h2osno [kg/m^2]
    real(rk8) :: h2osno_temp
    ! sum of snow ice if snow layers found above unfrozen lake [kg/m&2]
    real(rk8) :: sumsnowice(lbc:ubc)
    ! true if top lake layer is unfrozen with snow layers above
    logical  :: unfrozen(lbc:ubc)
    real(rk8) :: heatrem           ! used in case above [J/m^2]
    real(rk8) :: heatsum(lbc:ubc)  ! used in case above [J/m^2]
    real(rk8) :: snowmass          ! liquid+ice snow mass in a layer [kg/m2]
    ! temporary factor used to correct for snow capping
    real(rk8) :: snowcap_scl_fct
    ! assumed snow bulk density (for lakes w/out resolved snow layers) [kg/m^3]
    ! Should only be used for frost below.
    real(rk8), parameter :: snow_bd = 250._rk8

    ! Assign local pointers to derived subtypes components (gridcell-level)

    forc_rain => clm_a2l%forc_rain
    forc_snow => clm_a2l%forc_snow
    forc_t    => clm_a2l%forc_t

    ! Assign local pointers to derived subtypes components (column-level)

    frac_sno_eff  => clm3%g%l%c%cps%frac_sno_eff
    cgridcell     => clm3%g%l%c%gridcell
    clandunit     => clm3%g%l%c%landunit
    snl           => clm3%g%l%c%cps%snl
    t_grnd        => clm3%g%l%c%ces%t_grnd
    h2osno        => clm3%g%l%c%cws%h2osno
    snowice       => clm3%g%l%c%cws%snowice
    snowliq       => clm3%g%l%c%cws%snowliq
    zwt           => clm3%g%l%c%cws%zwt
    fcov          => clm3%g%l%c%cws%fcov
    fsat           => clm3%g%l%c%cws%fsat
    qcharge       => clm3%g%l%c%cws%qcharge
    qflx_top_soil => clm3%g%l%c%cwf%qflx_top_soil
    qflx_prec_grnd_col => clm3%g%l%c%cwf%pwf_a%qflx_prec_grnd
    qflx_evap_grnd_col => clm3%g%l%c%cwf%pwf_a%qflx_evap_grnd
    qflx_dew_grnd_col  => clm3%g%l%c%cwf%pwf_a%qflx_dew_grnd
    qflx_dew_snow_col  => clm3%g%l%c%cwf%pwf_a%qflx_dew_snow
    qflx_sub_snow_col  => clm3%g%l%c%cwf%pwf_a%qflx_sub_snow
    qflx_snwcp_ice_col => clm3%g%l%c%cwf%pwf_a%qflx_snwcp_ice
    watsat        => clm3%g%l%c%cps%watsat
    z             => clm3%g%l%c%cps%z
    dz            => clm3%g%l%c%cps%dz
    zi            => clm3%g%l%c%cps%zi
    t_soisno      => clm3%g%l%c%ces%t_soisno
    h2osoi_ice    => clm3%g%l%c%cws%h2osoi_ice
    h2osoi_liq    => clm3%g%l%c%cws%h2osoi_liq
    h2osoi_vol    => clm3%g%l%c%cws%h2osoi_vol
    qflx_drain    => clm3%g%l%c%cwf%qflx_drain
    qflx_surf     => clm3%g%l%c%cwf%qflx_surf
    qflx_infl     => clm3%g%l%c%cwf%qflx_infl
    qflx_qrgwl    => clm3%g%l%c%cwf%qflx_qrgwl
    qflx_runoff    => clm3%g%l%c%cwf%qflx_runoff
    qflx_irrig    =>  clm3%g%l%c%cwf%qflx_irrig
    endwb         => clm3%g%l%c%cwbal%endwb
    begwb         => clm3%g%l%c%cwbal%begwb
    dz_lake       => clm3%g%l%c%cps%dz_lake
    t_lake         => clm3%g%l%c%ces%t_lake
    lake_icefrac   => clm3%g%l%c%cws%lake_icefrac
    do_capsnow    => clm3%g%l%c%cps%do_capsnow
    snow_depth        => clm3%g%l%c%cps%snow_depth
    qflx_snow_grnd_col => clm3%g%l%c%cwf%pwf_a%qflx_snow_grnd
    frac_iceold   => clm3%g%l%c%cps%frac_iceold
    qflx_evap_tot_col => clm3%g%l%c%cwf%pwf_a%qflx_evap_tot
    soilalpha    => clm3%g%l%c%cws%soilalpha
    zwt          => clm3%g%l%c%cws%zwt
    rootr_column => clm3%g%l%c%cps%rootr_column
    qflx_rain_grnd_col => clm3%g%l%c%cwf%pwf_a%qflx_rain_grnd
    qflx_snomelt => clm3%g%l%c%cwf%qflx_snomelt
    eflx_snomelt => clm3%g%l%c%cef%eflx_snomelt
    ! Use column variables here.
    qflx_snwcp_ice_col => clm3%g%l%c%cwf%pwf_a%qflx_snwcp_ice
    qflx_snwcp_liq_col => clm3%g%l%c%cwf%pwf_a%qflx_snwcp_liq
    !SNICAR variables from Hydrology1
    snw_rds            => clm3%g%l%c%cps%snw_rds
    mss_bcpho          => clm3%g%l%c%cps%mss_bcpho
    mss_bcphi          => clm3%g%l%c%cps%mss_bcphi
    mss_bctot          => clm3%g%l%c%cps%mss_bctot
    mss_bc_col         => clm3%g%l%c%cps%mss_bc_col
    mss_bc_top         => clm3%g%l%c%cps%mss_bc_top
    mss_ocpho          => clm3%g%l%c%cps%mss_ocpho
    mss_ocphi          => clm3%g%l%c%cps%mss_ocphi
    mss_octot          => clm3%g%l%c%cps%mss_octot
    mss_oc_col         => clm3%g%l%c%cps%mss_oc_col
    mss_oc_top         => clm3%g%l%c%cps%mss_oc_top
    mss_dst1           => clm3%g%l%c%cps%mss_dst1
    mss_dst2           => clm3%g%l%c%cps%mss_dst2
    mss_dst3           => clm3%g%l%c%cps%mss_dst3
    mss_dst4           => clm3%g%l%c%cps%mss_dst4
    mss_dsttot         => clm3%g%l%c%cps%mss_dsttot
    mss_dst_col        => clm3%g%l%c%cps%mss_dst_col
    mss_dst_top        => clm3%g%l%c%cps%mss_dst_top
    ! Diagnostics
    snot_top          => clm3%g%l%c%cps%snot_top
    dTdz_top          => clm3%g%l%c%cps%dTdz_top
    snw_rds_top       => clm3%g%l%c%cps%snw_rds_top
    sno_liq_top       => clm3%g%l%c%cps%sno_liq_top
    h2osno_top        => clm3%g%l%c%cps%h2osno_top
    ! SNICAR variables from Hydrology2
    mss_cnc_bcphi     => clm3%g%l%c%cps%mss_cnc_bcphi
    mss_cnc_bcpho     => clm3%g%l%c%cps%mss_cnc_bcpho
    mss_cnc_ocphi     => clm3%g%l%c%cps%mss_cnc_ocphi
    mss_cnc_ocpho     => clm3%g%l%c%cps%mss_cnc_ocpho
    mss_cnc_dst1      => clm3%g%l%c%cps%mss_cnc_dst1
    mss_cnc_dst2      => clm3%g%l%c%cps%mss_cnc_dst2
    mss_cnc_dst3      => clm3%g%l%c%cps%mss_cnc_dst3
    mss_cnc_dst4      => clm3%g%l%c%cps%mss_cnc_dst4
    ! Flooding terms
    qflx_floodg       => clm_a2l%forc_flood
    qflx_floodc       => clm3%g%l%c%cwf%qflx_floodc
    frost_table       => clm3%g%l%c%cws%frost_table
    zwt_perched       => clm3%g%l%c%cws%zwt_perched
    qflx_drain_perched=> clm3%g%l%c%cwf%qflx_drain_perched
    qflx_h2osfc_surf  => clm3%g%l%c%cwf%qflx_h2osfc_surf
    qflx_snow_melt    => clm3%g%l%c%cwf%qflx_snow_melt
    qflx_rsub_sat     => clm3%g%l%c%cwf%qflx_rsub_sat
    qflx_top_soil     => clm3%g%l%c%cwf%qflx_top_soil
    qflx_sl_top_soil  => clm3%g%l%c%cwf%qflx_sl_top_soil



    ! Assign local pointers to derived type members (pft-level)

    pcolumn       => clm3%g%l%c%p%column
    pgridcell      => clm3%g%l%c%p%gridcell
    qflx_sub_snow  => clm3%g%l%c%p%pwf%qflx_sub_snow
    qflx_evap_grnd => clm3%g%l%c%p%pwf%qflx_evap_grnd
    qflx_dew_snow  => clm3%g%l%c%p%pwf%qflx_dew_snow
    qflx_dew_grnd  => clm3%g%l%c%p%pwf%qflx_dew_grnd
    qflx_prec_grnd     => clm3%g%l%c%p%pwf%qflx_prec_grnd
    qflx_snow_grnd_pft => clm3%g%l%c%p%pwf%qflx_snow_grnd
    qflx_rain_grnd     => clm3%g%l%c%p%pwf%qflx_rain_grnd
    qflx_evap_tot => clm3%g%l%c%p%pwf%qflx_evap_tot
    qflx_evap_soi  => clm3%g%l%c%p%pwf%qflx_evap_soi
    qflx_snwcp_ice => clm3%g%l%c%p%pwf%qflx_snwcp_ice
    qflx_snwcp_liq => clm3%g%l%c%p%pwf%qflx_snwcp_liq
    eflx_sh_tot    => clm3%g%l%c%p%pef%eflx_sh_tot
    eflx_sh_grnd   => clm3%g%l%c%p%pef%eflx_sh_grnd
    eflx_soil_grnd => clm3%g%l%c%p%pef%eflx_soil_grnd
    eflx_gnet      => clm3%g%l%c%p%pef%eflx_gnet
    eflx_grnd_lake => clm3%g%l%c%p%pef%eflx_grnd_lake


    ! Add soil water to water balance.
    do concurrent ( fc = 1:num_lakec )
      c = filter_lakec(fc)
      !$acc loop seq
      do j = 1, nlevgrnd
        begwb(c) = begwb(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
      end do
    end do

    ! Do precipitation onto ground, etc., from Hydrology1.

    do concurrent ( fp = 1:num_lakep )
      p = filter_lakep(fp)
      g = pgridcell(p)
      c = pcolumn(p)

      qflx_prec_grnd_snow(p) = forc_snow(g)
      qflx_prec_grnd_rain(p) = forc_rain(g)
      qflx_prec_grnd(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)

      if (do_capsnow(c)) then
        qflx_snwcp_ice(p) = qflx_prec_grnd_snow(p)
        qflx_snwcp_liq(p) = qflx_prec_grnd_rain(p)
        qflx_snow_grnd_pft(p) = 0._rk8
        qflx_rain_grnd(p) = 0._rk8
      else
        qflx_snwcp_ice(p) = 0._rk8
        qflx_snwcp_liq(p) = 0._rk8
        ! ice onto ground (mm/s)
        qflx_snow_grnd_pft(p) = qflx_prec_grnd_snow(p)
        ! liquid water onto ground (mm/s)
        qflx_rain_grnd(p)     = qflx_prec_grnd_rain(p)
      end if
      ! Assuming one PFT; needed for below
      qflx_snow_grnd_col(c) = qflx_snow_grnd_pft(p)
      qflx_rain_grnd_col(c) = qflx_rain_grnd(p)
    end do ! (end pft loop)

    ! Determine snow height and snow water

    do concurrent ( fc = 1:num_lakec )
      c = filter_lakec(fc)
      g = cgridcell(c)

      ! Use Alta relationship, Anderson(1976); LaChapelle(1961),
      ! U.S.Department of Agriculture Forest Service, Project F,
      ! Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.

      if (do_capsnow(c)) then
        dz_snowf = 0._rk8
      else
        if (forc_t(g) > tfrz + 2._rk8) then
          bifall=50._rk8 + 1.7_rk8*(17.0_rk8)**1.5_rk8
        else if (forc_t(g) > tfrz - 15._rk8) then
          bifall=50._rk8 + 1.7_rk8*(forc_t(g) - tfrz + 15._rk8)**1.5_rk8
        else
          bifall=50._rk8
        end if
        dz_snowf = qflx_snow_grnd_col(c)/bifall
        snow_depth(c) = snow_depth(c) + dz_snowf*dtsrf
        ! snow water equivalent (mm)
        h2osno(c) = h2osno(c) + qflx_snow_grnd_col(c)*dtsrf
      end if

      ! When the snow accumulation exceeds 40 mm, initialize snow layer
      ! Currently, the water temperature for the precipitation is simply set
      ! as the surface air temperature

      newnode = 0    ! flag for when snow node will be initialized
      if ( snl(c) == 0 .and. &
           qflx_snow_grnd_col(c) > 0.0_rk8 .and. &
           snow_depth(c) >= 0.01_rk8 + lsadz) then
        newnode = 1
        snl(c) = -1
        dz(c,0) = snow_depth(c)                       ! meter
        z(c,0) = -0.5_rk8*dz(c,0)
        zi(c,-1) = -dz(c,0)
        t_soisno(c,0) = min(tfrz, forc_t(g))      ! K
        h2osoi_ice(c,0) = h2osno(c)               ! kg/m2
        h2osoi_liq(c,0) = 0._rk8                   ! kg/m2
        frac_iceold(c,0) = 1._rk8

        ! intitialize SNICAR variables for fresh snow:
        snw_rds(c,0)    = snw_rds_min

        !$acc loop seq
        do j = -nlevsno+1, 0
          mss_bcpho(c,j)  = 0._rk8
          mss_bcphi(c,j)  = 0._rk8
          mss_bctot(c,j)  = 0._rk8
          mss_ocpho(c,j)  = 0._rk8
          mss_ocphi(c,j)  = 0._rk8
          mss_octot(c,j)  = 0._rk8
          mss_dst1(c,j)   = 0._rk8
          mss_dst2(c,j)   = 0._rk8
          mss_dst3(c,j)   = 0._rk8
          mss_dst4(c,j)   = 0._rk8
          mss_dsttot(c,j) = 0._rk8
        end do
        mss_bc_col(c)   = 0._rk8
        mss_bc_top(c)   = 0._rk8
        mss_oc_col(c)   = 0._rk8
        mss_oc_top(c)   = 0._rk8
        mss_dst_col(c)  = 0._rk8
        mss_dst_top(c)  = 0._rk8
      end if

      ! The change of ice partial density of surface node due to precipitation.
      ! Only ice part of snowfall is added here, the liquid part will be added
      ! later.

      if (snl(c) < 0 .and. newnode == 0) then
        h2osoi_ice(c,snl(c)+1) = h2osoi_ice(c,snl(c)+1) + &
                dtsrf*qflx_snow_grnd_col(c)
        dz(c,snl(c)+1) = dz(c,snl(c)+1)+dz_snowf*dtsrf
      end if

    end do

    ! Calculate sublimation and dew, adapted from HydrologyLake
    ! and Biogeophysics2.
    do concurrent ( fp = 1:num_lakep )
      qflx_evap_grnd(p) = 0._rk8
      qflx_sub_snow(p) = 0._rk8
      qflx_dew_snow(p) = 0._rk8
      qflx_dew_grnd(p) = 0._rk8
    end do

    do fp = 1,num_lakep
      p = filter_lakep(fp)
      c = pcolumn(p)
      jtop = snl(c)+1

      if (jtop <= 0) then ! snow layers
        j = jtop
        ! Assign ground evaporation to sublimation from soil ice or to dew
        ! on snow or ground

        if (qflx_evap_soi(p) >= 0._rk8) then
          ! for evaporation partitioning between liquid evap and ice
          ! sublimation, use the ratio of liquid to (liquid+ice) in the
          ! top layer to determine split Since we're not limiting evap
          ! over lakes, but still can't remove more from top
          ! snow layer than there is there, create temp. limited evap_soi.
          qflx_evap_soi_lim = min(qflx_evap_soi(p), &
                  (h2osoi_liq(c,j)+h2osoi_ice(c,j))/dtsrf)
          if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._rk8) then
            qflx_evap_grnd(p) = max(qflx_evap_soi_lim*(h2osoi_liq(c,j) / &
                    (h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._rk8)
          else
            qflx_evap_grnd(p) = 0._rk8
          end if
          qflx_sub_snow(p) = qflx_evap_soi_lim - qflx_evap_grnd(p)
        else
          ! if (t_grnd(c) < tfrz) then
          ! Causes rare blowup when thin snow layer should completely melt
          ! and has a high temp after thermal physics,
          ! but then is not eliminated in SnowHydrology because of this
          ! added frost. Also see below removal of
          ! completely melted single snow layer.
          if (t_grnd(c) < tfrz .and. t_soisno(c,j) < tfrz) then
            qflx_dew_snow(p) = abs(qflx_evap_soi(p))
            ! If top layer is only snow layer, SnowHydrology won't
            ! eliminate it if dew is added.
          else if (j < 0 .or. &
                  (t_grnd(c) == tfrz .and. t_soisno(c,j) == tfrz)) then
            qflx_dew_grnd(p) = abs(qflx_evap_soi(p))
          end if
        end if
        ! Update the pft-level qflx_snowcap
        ! This was moved in from Hydrology2 to keep all pft-level
        ! calculations out of Hydrology2
        if (do_capsnow(c)) then
          qflx_snwcp_ice(p) = qflx_snwcp_ice(p) + qflx_dew_snow(p)
          qflx_snwcp_liq(p) = qflx_snwcp_liq(p) + qflx_dew_grnd(p)
        end if

      else
        ! No snow layers: do as in HydrologyLake but with actual
        ! clmtype variables
        if (qflx_evap_soi(p) >= 0._rk8) then
          ! Sublimation: do not allow for more sublimation than there is snow
          ! after melt.  Remaining surface evaporation used for infiltration.
          qflx_sub_snow(p) = min(qflx_evap_soi(p), h2osno(c)/dtsrf)
          qflx_evap_grnd(p) = qflx_evap_soi(p) - qflx_sub_snow(p)
        else
          if (t_grnd(c) < tfrz-0.1_rk8) then
            qflx_dew_snow(p) = abs(qflx_evap_soi(p))
          else
            qflx_dew_grnd(p) = abs(qflx_evap_soi(p))
          end if
        end if

        ! Update snow pack for dew & sub.

        h2osno_temp = h2osno(c)
        if (do_capsnow(c)) then
          h2osno(c) = h2osno(c) - qflx_sub_snow(p)*dtsrf
          qflx_snwcp_ice(p) = qflx_snwcp_ice(p) + qflx_dew_snow(p)
          qflx_snwcp_liq(p) = qflx_snwcp_liq(p) + qflx_dew_grnd(p)
        else
          h2osno(c) = h2osno(c) + (-qflx_sub_snow(p)+qflx_dew_snow(p))*dtsrf
        end if
        if (h2osno_temp > 0._rk8) then
          snow_depth(c) = snow_depth(c) * h2osno(c) / h2osno_temp
        else
          !Assume a constant snow bulk density = 250.
          snow_depth(c) = h2osno(c)/snow_bd
        end if

#if (defined PERGRO)
        if (abs(h2osno(c)) < 1.e-10_rk8) h2osno(c) = 0._rk8
        if (h2osno(c) == 0._rk8) snow_depth(c) = 0._rk8
#else
        h2osno(c) = max(h2osno(c), 0._rk8)
#endif

      end if

      qflx_snwcp_ice_col(c) = qflx_snwcp_ice(p)
      qflx_snwcp_liq_col(c) = qflx_snwcp_liq(p)

    end do

    ! pft averages must be done here -- BEFORE SNOW CALCULATIONS AS THEY USE IT.
    ! for output to history tape and other uses
    ! (note that pft2col is called before SLakeHydrology, so we can't
    ! use that routine to do these column -> pft averages)
    do concurrent ( fp = 1:num_lakep )
      p = filter_lakep(fp)
      c = pcolumn(p)
      qflx_evap_tot_col(c) = qflx_evap_tot(p)
      qflx_prec_grnd_col(c) = qflx_prec_grnd(p)
      qflx_evap_grnd_col(c) = qflx_evap_grnd(p)
      qflx_dew_grnd_col(c) = qflx_dew_grnd(p)
      qflx_dew_snow_col(c) = qflx_dew_snow(p)
      qflx_sub_snow_col(c) = qflx_sub_snow(p)
    end do

    ! Determine initial snow/no-snow filters (will be modified possibly by
    ! routines CombineSnowLayers and DivideSnowLayers below)

    call BuildSnowFilter(lbc, ubc, num_lakec, filter_lakec, &
         num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, &
         filter_shlakenosnowc)

    ! specify snow fraction
    do concurrent ( fc = 1:num_lakec )
      c = filter_lakec(fc)
      if (h2osno(c) > 0.0_rk8) then
        frac_sno_eff(c)     = 1._rk8
      else
        frac_sno_eff(c)     = 0._rk8
      end if
    end do

    ! Determine the change of snow mass and the snow water onto soil

    call SnowWater(lbc, ubc, num_shlakesnowc, filter_shlakesnowc, &
            num_shlakenosnowc, filter_shlakenosnowc)

    ! Determine soil hydrology
    ! Here this consists only of making sure that soil is saturated even
    ! as it melts and pore space opens up. Conversely, if excess ice is
    ! melting and the liquid water exceeds the saturation value, then
    ! remove water.

    ! changed to nlevsoi on 8/11/10 to make consistent with non-lake bedrock
    do concurrent ( fc = 1:num_lakec, j = 1:nlevsoi )
      c = filter_lakec(fc)

      h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + &
                        h2osoi_ice(c,j)/(dz(c,j)*denice)
      ! Could have changed during phase change! (Added 8/11/10)

      if (h2osoi_vol(c,j) < watsat(c,j)) then
        h2osoi_liq(c,j) = (watsat(c,j)*dz(c,j) - &
                h2osoi_ice(c,j)/denice)*denh2o
        ! h2osoi_vol will be updated below, and this water addition
        ! will come from qflx_qrgwl
      else if (h2osoi_liq(c,j) > watsat(c,j)*denh2o*dz(c,j)) then
        h2osoi_liq(c,j) = watsat(c,j)*denh2o*dz(c,j)
        ! Another way to do this would be: if h2osoi_vol > watsat
        ! then remove min(h2osoi_liq, (h2osoi_vol-watsat)*dz*denh2o)
        ! from h2osoi_liq.  The question is whether the excess ice
        ! melts first or last (or simultaneously) to the pore ice.
        ! Because excess ice is often in chunks, requiring greater
        ! convergence of heat to melt, assume it melts last.
        ! This will also improve the initialization behavior or in
        ! an occasionally warm year, the excess ice
        ! won't start going away if a layer is briefly at freezing.

        ! Allow up to 10% excess ice over watsat in refreezing soil,
        ! e.g. heaving soil.  (As with > 10% excess ice modeling,
        ! and for the lake water, the thermal conductivity will be
        ! adjusted down to compensate for the fact that the nominal
        ! dz is smaller than the real soil volume.)
        !  The current solution is consistent but perhaps unrealistic
        ! in real soils, where slow drainage may occur during freezing;
        ! drainage is only assumed to occur here when >10% excess
        ! ice melts. The latter is more likely to be permanent rather
        ! than seasonal anyway. Attempting to remove the
        ! ice volume after some has already frozen during the timestep
        ! would not conserve energy unless this were
        ! incorporated into the ice stream.
      end if
    end do

    ! Natural compaction and metamorphosis.

    call SnowCompaction(lbc, ubc, num_shlakesnowc, filter_shlakesnowc)

    ! Combine thin snow elements

    call CombineSnowLayers(lbc, ubc, num_shlakesnowc, filter_shlakesnowc)

    ! Divide thick snow elements

    call DivideSnowLayers_Lake(lbc, ubc, num_shlakesnowc, filter_shlakesnowc)

    ! Check for single completely unfrozen snow layer over lake.
    ! Modeling this ponding is unnecessary and can cause instability
    ! after the timestep when melt is completed, as the temperature
    ! after melt can be excessive because the fluxes were calculated
    ! with a fixed ground temperature of freezing, but the
    ! phase change was unable to restore the temperature to freezing.
    do concurrent ( fp = 1:num_lakep )
      p = filter_lakep(fp)
      c = pcolumn(p)
      j = 0
      if (snl(c) == -1 .and. h2osoi_ice(c,j) == 0._rk8) then
        ! Remove layer
        ! Take extra heat of layer and release to sensible heat in
        ! order to maintain energy conservation.
        heatrem = cpliq*h2osoi_liq(c,j)*(t_soisno(c,j) - tfrz)
        eflx_sh_tot(p) = eflx_sh_tot(p) + heatrem/dtsrf
        ! Added this line 7/22/11 for consistency.
        eflx_sh_grnd(p) = eflx_sh_grnd(p) + heatrem/dtsrf
        eflx_soil_grnd(p) = eflx_soil_grnd(p) - heatrem/dtsrf
        eflx_gnet(p) = eflx_gnet(p) - heatrem/dtsrf
        eflx_grnd_lake(p) = eflx_grnd_lake(p) - heatrem/dtsrf
        qflx_sl_top_soil(c) = qflx_sl_top_soil(c) + h2osno(c)
        snl(c) = 0
        h2osno(c) = 0._rk8
        snow_depth(c) = 0._rk8
        ! Rest of snow layer book-keeping will be done below.
      end if
    end do

    ! Check for snow layers above lake with unfrozen top layer.  Mechanically,
    ! the snow will fall into the lake and melt or turn to ice.
    !  If the top layer has sufficient heat to melt the snow without
    ! freezing, then that will be done.
    ! Otherwise, the top layer will undergo freezing, but only if the
    ! top layer will not freeze completely.  Otherwise, let the snow
    ! layers persist and melt by diffusion.
    do concurrent ( fc = 1:num_lakec )
      c = filter_lakec(fc)
      if ( t_lake(c,1) > tfrz .and. &
           lake_icefrac(c,1) == 0._rk8 .and. snl(c) < 0) then
        unfrozen(c) = .true.
      else
        unfrozen(c) = .false.
      end if
    end do

    do concurrent ( fc = 1:num_lakec )
      c = filter_lakec(fc)
      !$acc loop seq
      do j = -nlevsno+1, 0
        if (unfrozen(c)) then
          if (j == -nlevsno+1) then
            sumsnowice(c) = 0._rk8
            heatsum(c) = 0._rk8
          end if
          if (j >= snl(c)+1) then
            sumsnowice(c) = sumsnowice(c) + h2osoi_ice(c,j)
            heatsum(c) = heatsum(c) + &
                    h2osoi_ice(c,j)*cpice*(tfrz - t_soisno(c,j)) + &
                    h2osoi_liq(c,j)*cpliq*(tfrz - t_soisno(c,j))
          end if
        end if
      end do
    end do

    do concurrent ( fc = 1:num_lakec )
      c = filter_lakec(fc)

      if (unfrozen(c)) then
        heatsum(c) = heatsum(c) + sumsnowice(c)*hfus
        heatrem = (t_lake(c,1) - tfrz)*cpliq*denh2o*dz_lake(c,1) - heatsum(c)

        if (heatrem + denh2o*dz_lake(c,1)*hfus > 0._rk8) then
          ! Remove snow and subtract the latent heat from the top layer.
          qflx_snomelt(c) = qflx_snomelt(c) + h2osno(c)/dtsrf
          eflx_snomelt(c) = eflx_snomelt(c) + h2osno(c)*hfus/dtsrf
          ! update snow melt for this case
          qflx_snow_melt(c)     = qflx_snow_melt(c)  + qflx_snomelt(c)

          qflx_sl_top_soil(c) = qflx_sl_top_soil(c) + h2osno(c)

          h2osno(c) = 0._rk8
          snow_depth(c) = 0._rk8
          snl(c) = 0
          ! The rest of the bookkeeping for the removed snow will be done below.
          if (heatrem > 0._rk8) then ! simply subtract the heat from the layer
            t_lake(c,1) = t_lake(c,1) - heatrem/(cpliq*denh2o*dz_lake(c,1))
          else !freeze part of the layer
            t_lake(c,1) = tfrz
            lake_icefrac(c,1) = -heatrem/(denh2o*dz_lake(c,1)*hfus)
          end if
        end if
      end if
    end do

    ! Set empty snow layers to zero

    do concurrent ( fc = 1:num_shlakesnowc, j = -nlevsno+1:0 )
      c = filter_shlakesnowc(fc)
      if (j <= snl(c) .and. snl(c) > -nlevsno) then
        h2osoi_ice(c,j) = 0._rk8
        h2osoi_liq(c,j) = 0._rk8
        t_soisno(c,j) = 0._rk8
        dz(c,j) = 0._rk8
        z(c,j) = 0._rk8
        zi(c,j-1) = 0._rk8
      end if
    end do

    ! Build new snow filter

    call BuildSnowFilter(lbc, ubc, num_lakec, filter_lakec, &
         num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, &
         filter_shlakenosnowc)

    ! Vertically average t_soisno and sum of h2osoi_liq and h2osoi_ice
    ! over all snow layers for history output

    do concurrent ( fc = 1:num_lakec )
      c = filter_lakec(fc)
      snowice(c) = 0._rk8
      snowliq(c) = 0._rk8
    end do

    do concurrent ( fc = 1:num_shlakesnowc )
      c = filter_shlakesnowc(fc)
      !$acc loop seq
      do j = -nlevsno+1, 0
        if (j >= snl(c)+1) then
          snowice(c) = snowice(c) + h2osoi_ice(c,j)
          snowliq(c) = snowliq(c) + h2osoi_liq(c,j)
        end if
      end do
    end do

    ! Determine ending water balance and volumetric soil water

    do concurrent ( fc = 1:num_lakec )
      c = filter_lakec(fc)
      endwb(c) = h2osno(c)
      !$acc loop seq
      do j = 1, nlevgrnd
        endwb(c) = endwb(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
        h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + &
                          h2osoi_ice(c,j)/(dz(c,j)*denice)
      end do
    end do

    ! Do history variables and set special landunit runoff (adapted
    ! from end of HydrologyLake)

    do concurrent ( fp = 1:num_lakep )
      p = filter_lakep(fp)
      c = pcolumn(p)
      g = pgridcell(p)

      zwt_perched(c)    = spval
      frost_table(c)    = spval
      qflx_drain_perched(c)= 0._rk8
      qflx_h2osfc_surf(c)  = 0._rk8
      qflx_rsub_sat(c)     = 0._rk8
      qflx_infl(c)      = 0._rk8
      qflx_surf(c)      = 0._rk8
      qflx_drain(c)     = 0._rk8
      qflx_irrig(c)     = 0._rk8
      soilalpha(c)      = spval
      zwt(c)            = spval
      fcov(c)           = spval
      fsat(c)           = spval
      qcharge(c)        = spval
      !$acc loop seq
      do j = 1, nlevgrnd
        rootr_column(c,j) = spval
      end do

      ! Insure water balance using qflx_qrgwl
      qflx_qrgwl(c)     = forc_rain(g) + forc_snow(g) - &
              qflx_evap_tot(p) - qflx_snwcp_ice(p) - &
              (endwb(c)-begwb(c))/dtsrf + qflx_floodg(g)
      qflx_floodc(c)    = qflx_floodg(g)
      qflx_runoff(c)    = qflx_drain(c) + qflx_surf(c) + qflx_qrgwl(c)
      qflx_top_soil(c)  = qflx_prec_grnd_rain(p) + qflx_snomelt(c)
    end do

    !  SNICAR Code and diagnostics

    !  Calculate column-integrated aerosol masses, and
    !  mass concentrations for radiative calculations and output
    !  (based on new snow level state, after SnowFilter is rebuilt.
    !  NEEDS TO BE AFTER SnowFiler is rebuilt, otherwise there
    !  can be zero snow layers but an active column in filter)

    do concurrent ( fc = 1:num_shlakesnowc )
      c = filter_shlakesnowc(fc)

      ! Zero column-integrated aerosol mass before summation
      mss_bc_col(c)  = 0._rk8
      mss_oc_col(c)  = 0._rk8
      mss_dst_col(c) = 0._rk8

      !$acc loop seq
      do j = -nlevsno+1, 0

        ! layer mass of snow:
        snowmass = h2osoi_ice(c,j)+h2osoi_liq(c,j)

        ! Correct the top layer aerosol mass to account for snow capping.
        ! This approach conserves the aerosol mass concentration
        ! (but not the aerosol amss) when snow-capping is invoked

        if (j == snl(c)+1) then
          if (do_capsnow(c)) then
            snowcap_scl_fct = snowmass / &
                (snowmass+(qflx_snwcp_ice_col(c)*dtsrf))
            ! Make sure column variable here
            mss_bcpho(c,j) = mss_bcpho(c,j)*snowcap_scl_fct
            mss_bcphi(c,j) = mss_bcphi(c,j)*snowcap_scl_fct
            mss_ocpho(c,j) = mss_ocpho(c,j)*snowcap_scl_fct
            mss_ocphi(c,j) = mss_ocphi(c,j)*snowcap_scl_fct

            mss_dst1(c,j)  = mss_dst1(c,j)*snowcap_scl_fct
            mss_dst2(c,j)  = mss_dst2(c,j)*snowcap_scl_fct
            mss_dst3(c,j)  = mss_dst3(c,j)*snowcap_scl_fct
            mss_dst4(c,j)  = mss_dst4(c,j)*snowcap_scl_fct
          end if
        end if

        if (j >= snl(c)+1) then
          mss_bctot(c,j)     = mss_bcpho(c,j) + mss_bcphi(c,j)
          mss_bc_col(c)      = mss_bc_col(c)  + mss_bctot(c,j)
          mss_cnc_bcphi(c,j) = mss_bcphi(c,j) / snowmass
          mss_cnc_bcpho(c,j) = mss_bcpho(c,j) / snowmass

          mss_octot(c,j)     = mss_ocpho(c,j) + mss_ocphi(c,j)
          mss_oc_col(c)      = mss_oc_col(c)  + mss_octot(c,j)
          mss_cnc_ocphi(c,j) = mss_ocphi(c,j) / snowmass
          mss_cnc_ocpho(c,j) = mss_ocpho(c,j) / snowmass

          mss_dsttot(c,j)    = mss_dst1(c,j)  + mss_dst2(c,j) + &
                  mss_dst3(c,j) + mss_dst4(c,j)
          mss_dst_col(c)     = mss_dst_col(c) + mss_dsttot(c,j)
          mss_cnc_dst1(c,j)  = mss_dst1(c,j)  / snowmass
          mss_cnc_dst2(c,j)  = mss_dst2(c,j)  / snowmass
          mss_cnc_dst3(c,j)  = mss_dst3(c,j)  / snowmass
          mss_cnc_dst4(c,j)  = mss_dst4(c,j)  / snowmass

        else
          !set variables of empty snow layers to zero
          snw_rds(c,j)       = 0._rk8

          mss_bcpho(c,j)     = 0._rk8
          mss_bcphi(c,j)     = 0._rk8
          mss_bctot(c,j)     = 0._rk8
          mss_cnc_bcphi(c,j) = 0._rk8
          mss_cnc_bcpho(c,j) = 0._rk8

          mss_ocpho(c,j)     = 0._rk8
          mss_ocphi(c,j)     = 0._rk8
          mss_octot(c,j)     = 0._rk8
          mss_cnc_ocphi(c,j) = 0._rk8
          mss_cnc_ocpho(c,j) = 0._rk8

          mss_dst1(c,j)      = 0._rk8
          mss_dst2(c,j)      = 0._rk8
          mss_dst3(c,j)      = 0._rk8
          mss_dst4(c,j)      = 0._rk8
          mss_dsttot(c,j)    = 0._rk8
          mss_cnc_dst1(c,j)  = 0._rk8
          mss_cnc_dst2(c,j)  = 0._rk8
          mss_cnc_dst3(c,j)  = 0._rk8
          mss_cnc_dst4(c,j)  = 0._rk8
        end if
      end do

      ! top-layer diagnostics
      h2osno_top(c)  = h2osoi_ice(c,snl(c)+1) + h2osoi_liq(c,snl(c)+1)
      mss_bc_top(c)  = mss_bctot(c,snl(c)+1)
      mss_oc_top(c)  = mss_octot(c,snl(c)+1)
      mss_dst_top(c) = mss_dsttot(c,snl(c)+1)
    end do

    ! Zero mass variables in columns without snow
    do concurrent ( fc = 1:num_shlakenosnowc )
      c = filter_shlakenosnowc(fc)

      h2osno_top(c)      = 0._rk8
      !$acc loop seq
      do j = -nlevsno+1, 0
        snw_rds(c,j)       = 0._rk8
        mss_bcpho(c,j)     = 0._rk8
        mss_bcphi(c,j)     = 0._rk8
        mss_bctot(c,j)     = 0._rk8
        mss_cnc_bcphi(c,j) = 0._rk8
        mss_cnc_bcpho(c,j) = 0._rk8
        mss_ocpho(c,j)     = 0._rk8
        mss_ocphi(c,j)     = 0._rk8
        mss_octot(c,j)     = 0._rk8
        mss_cnc_ocphi(c,j) = 0._rk8
        mss_cnc_ocpho(c,j) = 0._rk8
        mss_dst1(c,j)      = 0._rk8
        mss_dst2(c,j)      = 0._rk8
        mss_dst3(c,j)      = 0._rk8
        mss_dst4(c,j)      = 0._rk8
        mss_dsttot(c,j)    = 0._rk8
        mss_cnc_dst1(c,j)  = 0._rk8
        mss_cnc_dst2(c,j)  = 0._rk8
        mss_cnc_dst3(c,j)  = 0._rk8
        mss_cnc_dst4(c,j)  = 0._rk8
      end do

      mss_bc_top(c)      = 0._rk8
      mss_bc_col(c)      = 0._rk8
      mss_oc_top(c)      = 0._rk8
      mss_oc_col(c)      = 0._rk8
      mss_dst_top(c)     = 0._rk8
      mss_dst_col(c)     = 0._rk8

      ! top-layer diagnostics
      ! (spval is not averaged when computing history fields)
      snot_top(c)        = spval
      dTdz_top(c)        = spval
      snw_rds_top(c)     = spval
      sno_liq_top(c)     = spval
    end do

    ! Must be done here because the snow filter used in Hydrology2 &
    ! the Driver are for non-lake columns.
    call SnowAge_grain(lbc, ubc, num_shlakesnowc, filter_shlakesnowc, &
            num_shlakenosnowc, filter_shlakenosnowc)

  end subroutine SLakeHydrology

end module mod_clm_slakehydrology
! vim: tabstop=8 expandtab shiftwidth=2 softtabstop=2
